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Abstract. The shallow shelf approximation is a better "sliding law" for ice sheet mod- 
eling than those sliding laws in which basal velocity is a function of driving stress. The 
shallow shelf approximation as formulated by Schoof [2006a] is well-suited to this use. 
Our new thermomechanically coupled sliding scheme is based on a plasticity assump- 
tion about the strength of the saturated till underlying the ice sheet in which the till 
yield stress is given by a Mohr-Coulomb formula using a modeled pore water pressure. 
Using this scheme, our prognostic whole ice sheet model has convincing ice streams. Driv- 
ing stress is balanced in part by membrane stresses, the model is computable at high 
spatial resolution in parallel, it is stable with respect to parameter changes, and it pro- 
duces surface velocities seen in actual ice streams. 



1. Introduction 

A well-known difficulty with numerical ice sheet models 
is their inability to model the large range of ice flow speeds 
observed in real ice sheets [Shepard and Wingham, 2007; 
Truffer and Fahnestock, 2007; Vaughan and Arthern, 2007]. 
Observed surface speeds for ice flow in the Greenland ice 
sheet, for example, range from less than 10 meters per year 
in large areas of the interior [compare Greve, 1997b; Joughin 
et al, 1997] to more than 10 km per year in three outlet 
glaciers [Howat et al, 2007; Joughin et a/., 2004a]. Exist- 
ing Greenland ice sheet models have not, however, reported 
(published) ice surface speeds in excess of 100 m per year 
[Greve, 1997b, 2000; Ritz et al, 1997; Saito and Ahe-Ouchi, 
2005; Tarasov and Peltier, 2002]. 

Fast grounded ice flow, in ice streams and outlet glaciers 
to differing degrees [Truffer and Echelmeyer, 2003], arises 
from some combination of sliding, over a rigid or deformable 
mineral bed, and shear deformation of the lowest wet, dirty 
layers of ice. Unfortunately and fundamentally, however, re- 
mote sensing provides no high quality spatially-distributed 
observations of conditions at or near the ice base with which 
to constrain models of fast flow. There is a triple need to 
improve observations, to use existing surface observations 
more effectively, and to improve models of ice flow including 
sliding. 

This paper approaches modeling fast ice stream motion 
pragmatically, within the high-resolution, comprehensive, 
thermomechanically-coupled, and time-dependent Parallel 
Ice Sheet Model ["PISM"; Bueler et a/., 2008]. The basal 
mechanical model we add here is based on a spatially- 
distributed till friction angle [Paterson, 1994]. We demon- 
strate that our model responds in a reasonable way to 
changes in till friction angle and other major parameter 
choices including grid refinement. We believe that the model 
is a credible model of shallow ice streams. In our model, ice 
sheet geometry and thermodynamical fields (ice tempera- 
ture and effective thickness of basal water) evolve together 
within a unified shallow framework. 
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Fast-fiowing simulated ice is not useful in a model if it 
arises from unreasonable physics. All ice sheet models in- 
corporate approximations, and most models, including ours, 
use the actual shallowness of ice sheets to simplify the equa- 
tions and reduce computational cost. There are choices in 
parameterizing the sliding, however. We recall some con- 
tinuum flow models applied to ice sheets and sliding in the 
hierarchy in Figure 1. All the illustrated models describe ice 
as a slow, non-linear ly viscous, isotropic fluid, though these 
qualities are approximations too. 

The simplest and shallowest models are called the shal- 
low ice approximation (SIA) [Hutter, 1983; Morland and 
Johnson, 1980] and the shallow shelf approximation ["SSA"; 
Morland, 1987; Weis et al, 1999]. Rigorous small-parameter 
arguments explain how to simplify from Stokes to "higher- 
order models" [Blatter, 1995; Hindmarsh, 2004], and from 
the Blatter [1995] model to the SIA and SSA [Schoof and 
Hindmarsh, submitted]. 

Thermomechanically coupled, shallow, grounded, and 
non-sliding (frozen) base ice sheet models based upon the 
SIA are relatively well understood [Bueler et al, 2007; Payne 
et al, 2000]. Large portions of actual ice sheets have bases 
which experience minimal sliding and have modest bed to- 
pography. For those parts the nonsliding SIA is a second 
order [Fowler, 1997] theory which predicts a reasonable dis- 
tribution of flow at rates which compare well to observations 
(e.g. [Greve, 1997b]). 

As noted, faster ice flow is a combination of sliding over 
the mineral base along with deformation of an ice-and-till 
layer at the base. We necessarily lump these mechanisms as 
"sliding" in the language of this paper, because of the lack 
of observational techniques necessary to distinguish mecha- 
nisms at whole ice sheet scale. In any case, sliding applies 
a boundary force (stress) to the base of the ice mass, the 
effect of which is distributed by the stress balance. In our 
discussion of stress balance models below, components of 
the stress tensor will be separated into shear in planes par- 
allel to the geoid ( "horizontal plane shear" ) versus the other 
"membrane" stresses [Hindmarsh, 2004, 2006]. 

The SIA has no mechanism for balancing the driving 
stress partially by membrane stresses, nor any method of 
incorporating the sliding stress into a stress balance at all. 
Nonetheless "sliding laws" have been added to SIA models 
anyway [Greve, 1997b; Greve et al., 2006; Huybrechts and 
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de Wolde, 1999; Payne et al, 2000, among others]. Such 
laws necessarily describe the velocity of the base of the ice 
as a function of the driving stress, namely the product of the 
cryostatic pressure at the base times the surface slope [Pa- 
terson, 1994]. Furthermore, sliding is usually assumed to be 
insignificant when the ice base is frozen, but to "turn on" 
only when the ice base reaches the pressure-melting tem- 
perature [e.g. Greve et a/., 2006; Payne et al, 2000] or to 
change discontinuously at sufficient depth below sea level 
[Huybrechts and de Wolde, 1999; Tarasov and Peltier^ 2002]. 
Unfortunately, the underlying SIA continuum model then 
propagates jump discontinuities in the horizontal velocity 
field through the entire ice column. In turn, this produces 
unbounded vertical velocities because of the incompressibil- 
ity of the ice [Fowler^ 2001]. These facts about the con- 
tinuum model have the unfortunate consequence that nu- 
merical schemes for temperature or age will not converge 
under grid refinement, and, in fact, flow predictions from 
such numerical models become more unreasonable on finer 
grids (Appendix B). 

It remains critical, however, that modeled sliding depends 
on the ice temperature and/or on the amount of liquid water 
present at the ice base. Furthermore these quantities (tem- 
perature and basal water) must evolve if modeled ice streams 
are to exhibit the observed energy-balance-dependent be- 
havior actually seen [Raymond^ 2000; Schoof, 2004]. The 
coupling could potentially include a detailed model of till 
deformation [Pollard and DeConto, 2007], and perhaps also 
a distributed model for melt water conservation [Johnson 
and Fastook, 2002], but the model in this paper is simpler. 
An effective thickness for the water stored in the till is com- 
puted by time-integrating the rate of melt water production 
at the ice base. A negative melt rate (freeze-on) is allowed. 
The estimated pore water pressure is function of this ef- 
fective thickness. The pore water pressure is used when 
computing the effective pressure on the mineral till. The 
Mohr- Coulomb criterion then describes the yield stress of 
the saturated till as the product of the effective pressure on 
the till and the tangent of the till friction angle [Paters on, 
1994]. Till strength (yield stress) then enters as a term in 
a membrane stress balance, so the sliding velocity depends 
non-locally on till strength. 

Large ice shelves, with zero till strength, are well- 
described by the SSA model [MacAyeal et aL, 1996; Mor- 
land, 1987; Weis et a/., 1999]. Some published Antarctic 
ice sheet models use the SSA for the force balance in ice 
shelves [Huybrechts, 1990; Huybrechts and de Wolde, 1999; 
Ritz et al., 2001]. Furthermore, diagnostic models based on 
the "dragging ice shelf" extension of the SSA have been ap- 
plied to individual ice streams or ice stream basins [Hulbe 
and MacAyeal, 1999; MacAyeal, 1989]. These models have 
been exploited to recover the mechanical properties of the 
sliding ice base from observed surface velocities ["inverse 
modehng"; Joughin et al, 2001, 2004b]. Though MacAyeal 
[1989] proposed a linear viscous till deformation model, the 
mechanical properties of the till in real ice streams and 
outlet glaciers is a largely open question, with credible ar- 
guments for linear and plastic (Coulomb) models [Joughin 
et al., 2004b; Schoof , 2006b]. Power law formulations in 
between these extremes have also been considered [Schoof , 
submitted] . 

The current paper assumes that the till material behaves 
plasticly. This is because the best-understood continuum 
model which includes the SSA stress balance uses plastic till 
as part of a well-posed free boundary problem for the veloc- 
ity. In fact, Schoof [2006a] has recognized, in the isothermal 
and time- independent case, that including plastic failure of 
the basal till under flowing ice gives emergent ice streams 
within a whole ice sheet system. Our inclusion of Schoof 's 
model into a time- dependent model is new. The locations of 
sliding flow are not predetermined in Schoof's model, and 
ours. In ours these locations also evolve. 



We use the SSA model as a sliding law in a 
thermomechanically-coupled SIA model. The boundaries of 
the sliding regions are locations where the sliding velocities 
go to zero in a well-behaved way [Schoof , 2006a], so the 
problems described in Appendix B do not arise. The major- 
ity of the flow, by map-plane area or ice volume fraction, is 
by horizontal plane shear according to the non-sliding SIA. 
There are, however, ice streams which flow partly by the 
SIA but mostly by additional basal sliding constrained by 
the SSA balance of membrane stresses. The ice is allowed to 
slide anywhere, but sliding does not occur in the majority 
of the basal area because the till is too strong. Though till 
friction angle is time-independent, the actual till strength 
(yield stress) evolves in a thermomechanically coupled and 
time-dependent manner (through evolution of the modeled 
pore water pressure). 

The sliding velocity field arising at a particular time from 
the solution of the SSA is averaged with the velocity field 
which solves the nonsliding SIA at that time. This is the 
sense in which the SSA is a sliding law. As a result, the ve- 
locity field is a weighted average of two published ice sheet 
models, namely the nonsliding SIA [Hutter, 1983, etc.] and 
a plastic till form of the SSA [Schoof , 2006a]. The average 
is weighted toward the SIA if sliding is slow and toward the 
SSA if sliding is fast (subsection 2.8). 

As suggested by Figure 1, our stress balance combination 
model is not as complete as a number of "higher-order" or 
"full Stokes" (ISMIP-HOM) alternatives. The stress bal- 
ance for the Stokes model is, however, a three independent 
variable problem for the three components of the velocity 
field, at each time step. It has currently only been solved 
either for thermomechanically coupled but diagnostic prob- 
lems on glaciers [e.g. Zwinger et al., 2007] or with geometry 
evolution but isothermally and again on glacier-sized prob- 
lems [e.g. Pattyn et al., 2008]. The stress balance for the 
Blatter [1995] model is likewise a three independent vari- 
able problem at each time step which has been implemented 
on larger systems but not at high resolution [Pattyn, 2003; 
Saito et al, 2006]. Furthermore these higher order models 
have not demonstrated clear success for sliding flow even in 
flowline glacier modeling circumstances [Pattyn et al, 2008]. 

Our model is time-dependent in the usual sense that ice 
sheet geometry evolves according to a mass continuity equa- 
tion. The flow is thermomechanically coupled in the usual 
senses that the temperature affects the softness of the ice 
and there is dissipation heating from ice deformation and 
frictional heating from sliding. 

2. Continuum model 

Our description of the continuum model is organized into 
seven subsections: mass continuity, flow law, conservation 
of energy, basal melt, SIA version of the stress balance, SSA 
version of the stress balance, basal mechanics, and combina- 
tion of velocity fields. Notation follows Table 1. 

The time-independent boundary data are bed eleva- 
tion b(xi,X2), ice surface temperature Ts(xi, X2), accumu- 
lation rate M{xi,X2)- The modeled time-dependent un- 
knowns are the ice thickness H(t, xi,X2), surface elevation 
h{t, xi,X2), basal melt rate 5'(t,xi,X2), thickness of a melt 
water layer W{t, xi,X2), temperature T(t, xi,X2,X3), and 
vector velocity JJ{t,xi,X2,X3). The relation h = H -\- b 
holds always because only grounded ice sheets are con- 
sidered here. Let Dij be the strain rate tensor, that is, 
Dij = (1/2) {dUi/dxj ^dUj/dxi), aij the fufl (Cauchy) 
stress tensor, p = —(1/3) (an + (722 + (^33) the pressure, and 
Tij = CFij -h p5ij the deviatoric stress tensor. 

2.1. Mass continuity 

The upper surface z = h oi the ice is a free surface. Let 
Uh,i — Ui{z—h) be the surface value of the velocity. Then, 
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using a small slope assumption for the ice surface [Fowler, 
1997], 

dh 

-^=M + Uh,3-Vh- {Uh,uUh,2) (1) 

is the surface kinematical equation; V is in the horizontal 
variables. A similar equation applies at the base of the ice: 



= S + Ub,3-Vb-{Ub,i,Ub,2)- 



(2) 



Define Q = JJH, the horizontal ice fiux, where U is the 
vertically- averaged horizontal velocity. (For the non-sliding 
SIA (subsection 2.5 below), the expression Q = —Dsia V/i 
also applies, where Dsia is a positive diffusivity [Bueler 
et a/., 2007]. For membrane stress balanced sliding, a diffu- 
sivity can not be meaningfully defined because the ice fiux 
Q is not generally parallel to V/i. Instead we will treat 
the mass continuity problem associated to basal sliding as a 
generic transport problem (subsection 3.3 below).) 

Ice is an incompressible fiuid, so Dn -\- D22 + D33 = 0. 
The equation of mass continuity 



OH 



M-S-V 'Q. 



(3) 



can be derived from incompressibility by using the ice sur- 
face and base kinematic equations. We will solve Equation 
(3) numerically to compute the ice thickness at each grid 
point at each time step. 

The vertical velocity within the ice is, by incompressibil- 
ity, 



U3 



/X3 
V'{Ui,U2)dC 



(4) 



We have included the (ice-equivalent) basal melt rate S into 
the vertical velocity of the ice. In particular, the basal melt 
rate infiuences the vertical advection term in conservation 
of energy (subsection 2.3, next). Also, we have included the 
contribution of the basal melt rate S to the mass continuity 
equation (3); compare [Payne et al, 2000]. The basal melt 
rate will also feed back to affect ice dynamics by influencing 
the till yield stress (subsection 2.7). 

2.2. Flow law for ice 

Ice is modeled as a nonlinearly viscous isotropic fluid with 
a constitutive relation of Arrhenius- Glen- Nye form [Pater- 
son, 1994] 



(5) 



where T* = T -\- f3{h — z) is the pressure-corrected tem- 
perature, n = 3, and 2t^ = '^ij'Tij determines the second 
stress invariant r [Fowler, 1997]. This form of the flow law 
is used of the SIA stress balance (subsection 2.5). The ice 
softness (flow factor) A{T*) is determined by the formulas 
of Paterson and Budd [1982]: 



A(T*) 



(3.61 X 10-13) e-6 0xio^/(«^*), T* < 263.15K, 
(1.73 X 10+3) e-i3-9xioV(flT-)^ > 263.15K. 



For the SSA stress balance (subsection 2.6) we state the 
flow law in viscosity form 



Tij = 2v{T*,D)Dij. 



(6) 



In this case the nonlinear viscosity v satisfies 2i/(T*,D) = 
^Yiere 2D'^ = DijDij and 5(T*) = 
A(T*)-^/^ is the ice hardness. 

2.3. Conservation of energy 

We use a standard conservation of energy model for cold 
ice [not polythermal; Greve, 1997a]. Shallowness, namely a 



small ratio e = [^]/[^] of typical ice thickness [H] to typ- 
ical ice sheet width [/], simplifies the equation so that the 
horizontal conduction terms drop out. This simplification 
occurs both in the shallowness argument that leads to the 
SIA [Fowler, 1997] and the one that leads to the SSA [Weis 
et al, 1999], so we adopt it throughout: 



PiCi 



f + t/rf + + ^3^ 

Ot ox I OX2 OZ 



fe|^+E (7) 



Here pi,Ci, ki are constant material properties of the ice (see 
Table 1). The term in parentheses on the left of (7) is the 
material derivative dT/dt. We denote by E the rate at which 
straining is converted to heat and is applied to warming the 
ice volumetrically. 

In fact, recalling that ice sheets are a slow fiow, all work 
done by the driving stress of gravity is immediately de- 
posited as heat into the ice or is used for producing melt 
water at the base (subsection 2.7). In the case of the full 
Stokes model. 



: 2B{T*)D 



(l/n) + l 



(8) 



[Greve, 1997a]. We compute the strain heating S in equa- 
tion (7) using the shallow approximations of Eps described 
in subsection 2.8. 

Frictional heating occurs at the ice-lithosphere interface 
where the ice is sliding. The rate of heating is —Th • U^, a 
positive value with SI units of Jm~^ s~^ [e.g., chapter 10 of 
Paterson, 1994]. This rate times the area of the horizon- 
tal face of a grid cell times the length of a time step is an 
amount of heat which is added at each time step to the grid 
cell at the ice-lithosphere interface. 

Lithosphere "thermal inertia" is an important physical 
process in the thermomechanical regulation of ice stream 
fiow, and it is a standard part of paleo ice sheet modeling 
(e.g. [Greve, 1997b; Ritz et al, 1997]). We track the temper- 
ature of the lithosphere through a simplified conservation of 
energy model which again lacks horizontal conduction terms 
because of shallowness: 



q2ji 

K. — . (9) 



Here kr = 3.0 W(K m)-\ pr 



3300 kg m"^, and Cr = 
2009 J (kg K)"\ so Kr = kr/{prCr) 9.09 X 10"'^m^s"^ 
The thermal model extends downwards a distance of Dl = 
515 m into the lithosphere. This depth Dl is chosen for con- 
venience because no theory known to the authors identifies 
a preferred value. We observe that ice stream dynamics in 
our model are significantly different if Dl = 0, becayse with 
no lithosphere thermal model there is more oscillation in 
the streaming fiow, including more frequent complete halts 
of the fiow. 

The boundary conditions for Equation (7) apply a surface 
temperature to the top of each ice column and a geothermal 
fiux to the base of the lithosphere layer. 



T|(_,)=T.(xi,X2), 



2.4. Basal melt 



-k 



dT 



dx: 



3 \{z=b-DL) 



= Go. 



We compute a basal melt rate, and model the storage of 
melted water locally at the base of the ice column, as follows: 
If at some time step the ice in a grid cell reaches the pres- 
sure melting temperature Tq — Tq — f3{h — z) then our model 



X- 4 



BUELER AND BROWN: SHALLOW SHELF APPROXIMATION AS A "SLIDING LAW' 



converts a small fraction of the excess heat produced during 
that time step to melt water which appears at the base of 
the ice column. This fraction depends on the height above 
the bed. Concretely, for a time step At, "cold ice" Equation 
(7) predicts tentative temperature T{t + At), which might 
actually be above the pressure melting temperature. Our 
basal melt rate is 



nb+lOO ^ 

SAt=^ J (T{t + At)-nj 



6+100 



(10) 

The unitless number in square parentheses is the fraction of 
the melt water produced at height z which appears as free 
liquid water at the base. This fraction decreases linearly 
from 20% at the base to 0% at 100 meters above the base, 
and is zero above that. The basal melt rate ^S* can be neg- 
ative in the case where water is freezing on to the bottom 
of the ice column. Freeze-on only occurs if stored water is 
available at the base; storage is addressed next. 

The water produced by melting of ice near the base is 
locally stored in till. This water is described by an effective 
thickness W. It is available for weakening the till mechani- 
cally and for refreezing. Thickness W is updated by adding 
the basal melt rate ^S* at each time step, but with additional 
spatial diffusion: 



dW 

dt 



dxl 



dxl 



(11) 



where K^eit = L^/(2p), L = 20 km, and i = 1000 years. 
That is, the melt water is produced locally in each column 
according to (10) but it diffuses so that a delta mass spreads 
to a normal distribution with standard deviation of 20 km 
in 1000 years. Finally, W is limited to at least zero and 
at most Wo = 2 meters. The upper limit forms a minimal 
model for basal drainage. 

Once the basal melt rate ♦S' and the stored basal water 
thickness W are computed for time t + At, all temperatures 
within the ice are cut off at the pressure melting tempera- 
ture To*: T{t + At) = min{To*,f (t + At)}. 

A front-tracking [Greve, 1997a] or enthalpy gradient [As- 
chwanden and Blatter, submitted] polythermal model will 
improve the above melt-water production model. 

Subsection 2.7 addresses how W is used to compute till 
yield stress. The connection between stored basal water and 
till yield stress is a critical coupling in our model. 

2.5. SIA version of the stress balance 

The fundamental stress balance of glaciology is the Stokes 
model for a slow flowing fluid [Fowler, 1997]. Here "slow" 
has the precise meaning that inertia can be neglected in the 
force (stress) balance. In this paper we combine two shal- 
low versions of the Stokes model. One is the shallow ice 
approximation (SIA) [Hutter, 1983; Morland and Johnson, 
1980]. We use only the well-justified non-sliding version of 
the SIA, a non-Newtonian lubrication approximation for ice 
[Fowler, 1997]. The nonsliding SIA gives an ice sheet flow 
in which the driving force of gravity is balanced exclusively 
by horizontal plane shear. 

Denoting the velocity field which solves this stress bal- 
ance by u = {ui,U2,U3), the SIA can be written as a single 
equation 



In our model the integral in square brackets is computed 
numerically by the trapezoid rule on a vertical grid which is 
finest near the base of the ice. 

Throughout this paper the pressure is modeled as cyro- 
static, p = pg{h — z). This approximation holds in both the 
SIA [Fowler, 1997] and the SSA [Weis et al, 1999]. 

2.6. SSA version of the stress balance 

When there is significant sliding at the base the driving 
stress is significantly balanced by membrane stresses [Hind- 
marsh, 2006]. "Membrane stresses", called "longitudinal" 
in flowline models, are deviatoric components rii,r22,Ti2. 
Any membrane stress balance is nonlocal because the driv- 
ing stress is balanced by connection to distant ice either in 
the along-fiow direction or at the margins of regions of fast 
flow. The driving stress must be fully balanced by mem- 
brane stresses in cases where there is negligible traction at 
the base of the ice as in an ice shelf [Morland, 1987]. 

The simplest form of a membrane stress balance derivable 
from the Stokes model is the shallow shelf approximation 
["SSA" Weis et al., 1999]. It describes a fluid membrane of 
variable thickness in which gravity causes spreading. When 
applied to ice streams the membrane also experiences trac- 
tion (shear) forces on its base [MacAyeal, 1989]. 

We denote the velocity field which solves the SSA stress 
balance by v = {vi,V2,V3). The horizontal velocity {vi,V2) 
and the strain rates Dn, D22: D12 modeled by the SSA are 
independent of depth. Denote vertically- averaged hardness 
hy B, so BH = B^T"") dz. The vertically-averaged vis- 
cosity is 



B 



Dli + DI2 + D11D22 + J {D12 + D2if 



(14) 

[MacAyeal, 1989]. Following [Morland, 1987; Schoof, 2006a] 
we also define a vertically-integrated stress tensor 



T^j =2iyH (A, + {Dii + D22)Si^ 



(15) 



for i,j = 1,2. Note Tn + T22 is not generally zero so the 
membrane flow is conceptually "compressible" in the map- 
plane. In fact Til +T22 = (2zyif)(— 3D33) so a negative value 
for Til + T22 is a map-plane "pressure" causing an upward 
strain rate D33 = dvs/dz > 0. 

Let (r?,,i,r&,2) be the basal shear stress components ap- 
plied to the base of the ice (modeled in next subsection). 
The SSA is the pair of equations 



dTa ^ dTi2 ^ rr 

ox I 0x2 



dh 

dxi 



(16) 



for i = 1,2 [Schoof, 2006a]. A slightly more concrete form 
is [MacAyeal, 1989; Weis et al, 1999] 



d 

dxi 



29H 



^ dvi _^ dv2 
dxi dx2 



d 
dx2 



9H 



d 

dxi 



dvi ^ dv2 
dx2 dxi 



+ 



dx2 



2uH 



dvi ^ dv2 
dx2 dxi 



dvi _^ ^ dv2 
dxi dx2 



(17) 



dui du2 
dz ' dz 



-2{pgrA{T^){h- 



\Vhy-^Vh. (12) +Tb,2=pgH 



dh 
dx2 ' 



Equation (12) computes shear strain rates T>i3, D23 from lo- 
cal information, in the map-plane sense. Vertical integration 
gives 



[Ul,U2) 



-2{pgr iv/ir 



£a(t*)(/^- 



V/i. 
(13) 



Because our emphasis is on a novel interpretation of the 
SSA as a "sliding law", we note the heuristic which says 
Equation (17) is a splitting of the driving stress: 



(membrane stressesX 
held by viscous + 
deformation J 



^ stress held \ /a • • \ 
at base by = ' ' 
^till strength J 



\ stress J 
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In the extreme case of ice shelves the basal stress is zero 
so the left-most term fully supports the driving stress. In 
the other (singular) extreme case of strong till none of the 
driving stress is held by the membrane strain rates, in which 
case the SSA predicts no flow [Schoof, 2006a]. In that case 
we propose that the SIA should "take over" and predict flow 
by shear in horizontal planes. Again, this is the sense of our 
title. 

We will solve Equations (17) numerically at each time 
step to determine velocities Vi from evolving geometry and 
temperature. These equations are elliptic. More precisely, 
these equations can be derived from a variational principle 
applied to a convex and bounded-below functional [Equa- 
tion (3.13) of Schoof, 2006a], and the problem is well-posed. 
On the other hand, we will add dynamics, that is we will 
couple the SSA with Equations (3) and (7). Such an evolv- 
ing coupled system is not yet known to be mathematically 
well-posed. 

2.7. Basal mechanics 

We follow Schoof [2006a] and assume that the basal shear 
stress has plastic form 



Tb,i 



(^2 _^ ^2)1/2 



(18) 



for z = 1, 2, where Tc is a positive scalar yield stress [Pater- 
son, 1994]. Conceptually, because the till is assumed to be 
plastic, it supports applied stress without deformation un- 
til the applied stress equals the yield stress, at which point 
some amount of deformation occurs. 

In fact rb,i is only literally given by Equation (18) in map- 
plane locations where sliding is occurring. In locations where 
no sliding occurs, the driving stress on the right side of (17) 
is equal to the left side in the sense that the basal stress 
equals the driving stress. Where the base is not sliding, no 
division by zero occurs in the mathematically-precise form of 
Equations (17) because the problem is posed as minimizing 
a functional involving no division, though at its minimum 
the functional fails to have a well-defined gradient [Schoof, 
2006a, Equation (3.13)]. 

In our thermomechanically-coupled model, the value of 
the till yield stress is partly determined by the availabil- 
ity of basal water. The Mohr-Coulomb model for the yield 
stress of saturated till, which we adopt, is 



Tc = Co + (tan0) (pgH - p^j). 



(19) 



[Paterson, 1994, Chapter 8]. Here cq is the till cohesion, 
Pw is the pore water pressure (estimated below), and is a 
"till friction angle," a strength parameter for the till com- 
parable to "angle of repose" for granular piles. In this paper 
we take co = partly for simplicity and partly because Pa- 
terson [1994] notes an inconclusively broad observed range 
0-40 kPa for cq. The factor pgH—p^ is the effective pressure 
of the overlying ice on the mineral portion of the saturated 
till. The till is weakened by the presence of pressurized liq- 
uid water. 

In our model Tc necessarily represents the yield stress of 
the aggregate material at the base of an ice sheet, a poorly 
observed mixture of liquid water, ice, and granular till. 

Recall we model an effective thickness W of stored liquid 
water at the base of the ice column. This liquid water is in 
fact mixed with the solid part of the till, in some manner 
which we make no attempt to model. Nonetheless we use 
W to estimate the pore water pressure: 



p^ = 0.95 (pgH) 



W 
Wo 



(20) 



Note < W < 2 are imposed limits. We use Wo = 2 m 
throughout. Thus we model the pore water pressure locally 
as at most a fixed fraction (95%) of the ice overburden pres- 
sure pgH. When the base is frozen we have W = so pw = 



and the till is strong. With a till friction angle (f) = 15° and 
fully-saturated till, the yield stress Tc is 1.3% of overburden, 
while with = 2° the yield stress is 0.17% of overburden. 

2.8. Combining the velocities 

The use of an SSA result as the sliding velocity in an 
SIA model, the most novel aspect of our model, means com- 
bining the velocities computed by SIA and SSA to get the 
velocity used in mass continuity and conservation of energy. 
Recall that u = (1^1,1^2) is the SIA velocity (subsection 2.5) 
and V = {vi,V2) is the SSA velocity (subsection 2.6). We 
compute the combined horizontal velocity U = (C/i, C/2) by 



((7i,C/2) = /(|v|)u+(l-/(|v|))v, 
with |v|^ = Vi -\-vl and 

/(|v|) = l-^arctan(i^). 



(21) 



(22) 



The weighting function / has values between zero and one. 
See Figure 2. It satisfies the general requirements of smooth- 
ness, monotone decrease, /(|v|) ~ 1 for |v| small, and 
/(|v|) ~ for |v| large (e.g. |v| > 100 m/a). Many other 
choices satisfying these requirements are possible, and com- 
parison of trusted Stokes model results to those from this 
paper is a potential method for choosing an optimal weight- 
ing /(|v|). 

It follows from Equations (21) and (13) that 



Q = -/(|v|)DsiaV/i + (1 - /(|v|))vff, 



(23) 



where 



DsiA = 2{pgy 



nh 

\vh\^-^ J 



AiT'') {h - z)""^^ dz (24) 



is the positive diffusivity associated to the thermomechan- 
ical SIA [Bueler et al., 2007, and references therein]. Note 
that the flux is a linear combination of a vector pointing 
down the surface slope and a vector determined through 
the SSA stress balance. Only the former part of the flux 
is automatically diffusive for mass continuity Equation (3); 
compare [Pattyn, 2003]. We address the numerical approx- 
imation of Q and a time stepping scheme for Equation (3) 
in subsection 3.3. 

Equation (8) gives the strain heating rate Sps applicable 
to a full Stokes model. We approximate Eps by a com- 
putable and appropriately shallow quantity E. Strain rates 
for u and v do not sum in a simple way to yield the com- 
bined strain rates Dij for U, however. Nonetheless we pro- 
pose an approximate method for computing E in pieces, as 
follows. Denoting D{u)ij = (1/2) {dui/dxj -\- duj /dxi) and 
D{y)ij = (1/2) {^v^/^XJ + dvj/dxi), 

D^J ^ /(|v|)D(u)., + (1 - /(|v|)) i^(v).,. (25) 

The SIA says that each of D(u)ii, D(u)i2, i:>(u)2i, i:)(u)22 
are negligible D(v)i3, i^(v)3i, i^(v)23, ^(v)32 are all negli- 
gible in the SSA. Thus equation (25) can be rewritten 



(A,)^/(|v|) 



^0 D{u)i3^ 

D{U)23 

< J 



+ (l-/(|v|)) 



^ 


^(v)33> 
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Only entries on and above the diagonal need be displayed 
because of symmetry. It follows from (25) that 

= f{\v\f {D{u)ls + D{u)ls) + (1 - /(|v|))^ D{vf 

where Divf = D{w)l^ + D{w% + D{w)l2 + D(w)iiD{w)22. 
This yields a computable approximation of Sfs = 
25(r*)D(^+^)/^, in terms of quantities available in the SIA 
and SSA stress balance computations. 



3. Numerics and verification 

Each part of the continuum model of section 2 
is numerically approximated in PISM, an open source 
project (www.pism-docs.org). For the time-dependent, 
thermo-mechanically- coupled, non-sliding SIA, the numer- 
ical schemes used in PISM are described in [Bueler et a/., 
2007, Appendix A] . The PISM schemes described in this sec- 
tion are those additional schemes which were used to pro- 
duce results for this paper. 

All PISM numerical schemes share these properties: 

(i) they use finite difference approximations on a regular, 
rectangular grid in horizontal variables xi,X2, 

(ii) they do not use a rescaling of the vertical axis [e.g. as 
in Jenssen, 1977], but instead model the ice surface as a 
boundary between ice and air within the computational grid, 

(Hi) time stepping is adaptive because the stability con- 
ditions appropriate to each explicit or semi-implicit time- 
stepping numerical scheme are combined to give a next time- 
step based on the current geometry, velocity, and tempera- 
ture fields, and 

(iv) the PETSc library [Portable, Extensible Toolkit for 
Scientific computation; Balay and eight others^ 2006] man- 
ages parallel communication and the parallel, iterative solu- 
tion of linear systems. 

3.1. Approximation of the SSA 

Equations (17) are nonlinear in the velocities Vi. There- 
fore we solve them through an iteration in which the values 
of the velocities from the previous time step are used to com- 
pute an initial estimate of the vertically-integrated effective 
viscosity v^^^ The iteration then solves the linearized form 
of Equation (17) to compute new velocities and updated 
viscosity V^^^ with each iteration. The final values v from 
this iteration are the "SSA velocities" . Such an iteration has 
been widely used to solve the diagnostic problem of comput- 
ing the flow velocity in ice shelves [MacAyeal et al, 1996]. 
Appendix A details the PISM parallel implementation. 

We regularize Equations (14) and (18) to avoid divisions 
by zero. In fact, we must regularize the problem to even 
think of it as a system of partial differential equations, be- 
cause otherwise the correct formulation of the problem is a 
weaker form (a "variational inequality"). Following [Schoof, 
2006a] we introduce both an ice viscosity regularization and 
a till plasticity regularization. For the former, we replace 
Equation (14) by 



= ? [72 + ^11 + ^22 + D11D22 + ] {D12 + D2lf 
A 4 

-(26) 

where e has units of velocity and L^^ is a characteristic 
horizontal length scale for the membrane stresses. Our 
values t — \ meter per year and — 1000 km give 
e/Li, 3 X 10~^^s~^, a strain rate which should be com- 
pared to the typical membrane strain rates in a ice stream 
or shelf. These typical rates are of order 100 meters per 
year per 100 km, about 3 x 10~^^ s~^, in the Siple Coast ice 



streams \MacAyeal^ 1989] or the Ross ice shelf [MacAyeal 
et aL, 1996]. Because of the squares in Equation (26), for 
large parts of the sliding domain there is a 10^ difference in 
magnitudes between the regularizing constant (e/Lj^)^ and 
the largest of the other quantities in square brackets. 

For the plastic regularization we modify Equation (18) 
by choosing a small velocity S = 0.01 meters per year and 
defining ^ 

= (52+^2^^2)1/2 • (27) 

This regularization is only significant at or near locations 
where there is no sliding in the perfectly plastic model. De- 
scribed another way, we have made the till linearly viscous 
rb,i = —f^Vi^ with a very large coefficient [5 = rcS~^ , when 
sliding velocities are of order S or smaller. Thus the ice 
experiences slight sliding even where the till would not fail 
in the perfectly plastic limit. It follows from Equation (19) 
that this sliding is of order 6 where the base is melted but 
not yielding, and about 10~^^ where the base is frozen and 
not yielding. 

In the context of a higher-order Blatter [1995] model, 
Schoof [submitted] has recently shown that the ice flow 
velocity which results from the plasticity-regularized equa- 
tions converges to the purely-plastic result. We expect 
that a similar result applies in our vertically-integrated, 
thermomechanically-coupled, and geometry-evolving case. 

3.2. Verification of the SSA numerics 

We use an existing exact solution for primary verification 
of our numerical approximation of the "diagnostic" isother- 
mal SSA (including plastic till but without time evolution of 
stream geometry). "Verification" is used here as in compu- 
tational fluid dynamics [Wesseling , 2001] to mean the com- 
parison of numerical results to exact predictions of the con- 
tinuum model. 

The exact solution we use appears in section 4 of [Schoof, 
2006a], but it essentially arises in the analysis in [Raymond, 
2000] as well. Parameters and notation for the exact solu- 
tion are in Table 2. Because of translational symmetry in 
the downstream direction, the SSA with plastic till reduces 
to a one variable problem in the cross-stream coordinate y. 
The solution is a formula v{y) for the velocity which solves 
a free boundary problem describing a single ice stream with 
given till yield stress values Tx{y)- Because this is a free 
boundary problem, the width of the sliding ice stream is 
found simultaneously with the velocity profile v{y). Note 
that the thickness, bed slope, and surface slope are all con- 
stant. We use the particular instance illustrated in Figure 
3; it is called "Test I" in PISM [Bueler et al, 2008]. 

Our goal is to measure how well the finite difference 
scheme, described in Appendix A, solves this free bound- 
ary problem. Using a refinement path with Ay decreasing 
from 5 km to 40 m, the maximum and average numerical er- 
rors in velocity decay as shown in Figure 4. The numerical 
error at grid point k is \vnuni{yk) - '^exact (l/fc ) | • The max- 
imum of these pointwise errors converges rapidly down to 
a level of less than 10~^ m/a as the grid is refined. For 
comparison, the exact maximum velocity is 777.5 m/a. The 
average error is about three times smaller than the maxi- 
mum error. As shown in Figure 4, for "rough" grids with 
160 m < Ay < 5 km the errors decay at almost the optimal 
0{Ay^) rate for such finite difference schemes [Morton and 
Mayers, 2005] . For the finest two grids in Figure 4 the errors 
are significantly infiuenced by the viscosity iteration relative 
tolerance. This tolerance has been set to 5 x 10~^ and the 
linear iteration (Krylov solver) relative tolerance to 10~^^ 
for these verification experiments. 

With regard to this SSA free boundary problem, Schoof 
[2006a] also demonstrates that for any aspect ratio that 
could conceivably be regarded as "shallow", such as depth- 
to- width ratios of 0.1 or smaller, the velocity predictions of 



BUELER AND BROWN: SHALLOW SHELF APPROXIMATION AS A "SLIDING LAW' 



X- 7 



the SSA are very close to those of the full Stokes equations 
(applied to the same boundary values). 

3.3. Approximation of mass continuity 

The mass continuity Equation (3) is only strictly diffusive 
if the vertically-integrated flux Q is exactly in the downslope 
direction at all times and all locations on the ice sheet. But 
there is no reason to expect this to be true for real ice sheets. 
For instance, in ice streams with minimal bed topography 
and low surface slope it is likely that the flow direction is 
controlled largely by side drag and other effects unrelated to 
local surface slope. In our model some flux acts diffusively 
and some is treated as generic transport. 

In fact, recalling Equation (23), it follows that Equation 
(3) can be written 

^ = V • (^DVh?j - V • Vi7 - (V • v)H + (M - S). (28) 

where D — /(|v|)L)sia is a positive scalar and v = (1 — 
fiWl))^- Thus, temporarily ignoring the distinction be- 
tween H and /i. Equation (28) has diffusion, advection, re- 
action, and source terms on the right side (in that order). 

Our numerical approach to Equation (28) is to use the 
simplest explicit finite difference method which is condi- 
tionally stable and consistent [Morton and Mayers, 2005]. 
Though we describe the scheme next as though there were 
only one horizontal variable x, recovery of the (xi, X2) form 
is straightforward. We will use m in the superscript to de- 
note time step and j in the subscript to denote horizon- 
tal grid point xj . Also suppressing the superscript m on the 
right because all terms are evaluted at (i.e. explicitly), 
the scheme for (28) is 

- Hp ^ Dj+1/2 {hj+i - hj) - Dj-1/2 {hj - hj-i) 
At ~ Ax2 

(29) 

The weighted diffusivity i)j+i/2 is computed by the Mahaffy 
[1976] scheme as in [Bueler et al, 2007]. Our upwind scheme 
uses notation Up (/f.I'Dj) for Vj{Hj — Hj-i) if vj > and 
- Hj) if Vj < 0. 

Scheme (29) is conditionally stable. A condition applies 
to the explicit diffusion scheme: (max^ 1/2) < 0.5Ax^ 
[Morton and Mayers, 2005, chapter 2]. The CFL condition 
for the upwind scheme applies: (max^ vj^) At < Ax [Mor- 
ton and Mayers, 2005, chapter 4]. We adaptively take At to 
be the largest time step which satisfies both of these condi- 
tions at every point on the grid. Though the mass continu- 
ity equation is strongly nonlinear, stability follows because 
the maximum principle applies [Morton and Mayers, 2005, 
pages 16-18, 36-38, 50-51, and 94-96]. 

There is an additional condition for the stability of the 
conservation of energy scheme [Bueler et al, 2007], and it is 
also included in the adaptive time-stepping scheme. 

4. Experiments 

The model described above is not applied to a real ice 
sheet in this paper, but it is being applied to the whole 
Antarctic ice sheet in separate work (in preparation). Here 
we examine how the modeled flow speed and ice sheet vol- 
ume depend on certain parameters. We use four experiments 
which explore part of the parameter space and we study grid 
refinement. These experiments demonstrate the stability of 
the model with respect to parameter changes. 

Each of the four experiments starts from the same cir- 
cular ice sheet which is the initial state. It is very similar 



to the steady state (i.e. at 200 ka model years) of exper- 
iment A in EISMINT II [Payne et al, 2000], which is de- 
noted "EISIIA". The extent of the model domain and the 
net surface mass balance and surface temperature maps are 
all exactly as in EISIIA, and our initial state comes from a 
model run using only the non-sliding, thermomechanically- 
coupled SIA, as in EISIIA. Two differences are that ours 
includes a bedrock thermal model (Equation (9)) and mod- 
eled basal water effective thickness (Equation (11)). For our 
initial state, the basal melt rate is computed by Equation 
(10), and it is included in the mass continuity Equation (3) 
and in the vertical velocity (Equation (4)). Our initial state 
is on a 10 km horizontal grid, versus 25 km for EISIIA. 

The initial state and the experiments which follow use an 
unequally-spaced vertical grid in the ice (and air above, up 
to elevation 5000 m), with spacing increasing from 12.9 m 
at the bed to 87.1 m at the top. The vertical grid in bedrock 
has forty equal 12.9 m grid spaces. 

The numerically-computed initial state has radius about 
575 km, center thickness 3708.75 m, center (absolute) tem- 
perature at the ice base 256.24 K, and volume 2.208 x 
10^ km^. Figure 5 shows the thickness and the extent of 
melted base for the initial state. 

This initial ice sheet does not come from a sliding model, 
but it is "ready to slide" in locations where the till yield 
stress Tc is sufficiently low. On the one hand, if the till fric- 
tion angle were set to a large value everywhere (e.g. 15° 
or greater), then our model reduces to the non-sliding 
thermomechanically-coupled SIA with no sliding. On the 
other hand, if the till friction angle were set everywhere to 
a sufficiently weak value (e.g. 5°) then the sheet would slide 
in all radial directions. In fact, however, each of our four ex- 
periments includes a map of till friction angle with weak 
strips. Under most of the sheet (f) — 15°, while in the weak 
strips we set = 5°, except in experiment P4 which ex- 
plores this weak till parameter. The experiments start from 
the initial sheet described in the previous paragraphs, and 
run for 5000 model years. One run, experiment PI on a 10 
km grid, is extended to 100k model years. See Table 3. 

The spatial distribution of weak till for three of the four 
experiments is shown in Figure 6. Experiment PI has four 
weak till strips with widths 30, 50, 70, and 100 km, oriented 
in the cardinal directions. The cross stream profiles of till 
friction angle for the four strips are shown in Figure 7. 
Experiment P2 has three weak till strips each with width 70 
km, but oriented at 0°, 10°, and 45° relative to the closest 
cardinal grid direction. The goal of this experiment is to 
determine to what degree the grid alignment of the weak 
strip affects model outcome. Experiments P3 and P4 each 
have four weak strips in the cardinal directions with iden- 
tical widths 70 km. Experiment P3 has constant weakness 
of 5°, with cross-stream profile exactly as shown in the 70 
km width case in Figure 7, but there is differing bed slope 
in the four streams. We essentially use the "trough" bed 
topography from EISMINT II experiment I [Payne, 1997], 
but with total elevation drops of 0, 500, 1000, and 2000 m in 
the 650 km length of the troughs, giving slopes of 0, 0.077, 
0.154, and 0.308 percent, respectively. Figure 8 shows the 
map of bed elevation for P3. Finally, experiment P4 has dif- 
ferent values for the downstream value of (j). "Downstream" 
starts 400 km from the beginning of the weak strip. The dif- 
ferent strips have downstream values 2°, 3°, 8°, and 10°, 
respectively, while the upstream value is 5° for all strips. 

All experiments are grounded ice sheets without calv- 
ing fronts. The combined SIA-SSA model in this pa- 
per is applied at all points and there is no calving-front 
boundary condition. We have not analyzed the asymptotic 
shape of the resulting grounded margin for the combined 
model. The margin shape arising from the isothermal and 
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thermomechanically-coupled nonsliding SIA models is ana- 
lyzed in [Bueler et aL, 2005, 2007], while numerical scheme 
improvements for SIA margins appear in [Saito et a/., 2007], 
but these results do not apply directly to our case. 

In preparation for the next section, and to illustrate why 
the initial state is "ready to slide" in the weak till strips, 
in Figure 9 we compare the magnitude of the basal driving 
stress pgH\\/h\, for the initial state, to the till yield stress 
computed from Equation (19) for the till friction angle map 
of experiment PI. We see that the driving stress exceeds 
the yield stress only in the portions of the strips where the 
base has stored melt water; compare Figure 5. The frozen 
parts of the weak till strips, though they have till friction 
angle = 5°, actually have large yield stresses as a con- 
sequence of Equation (19). In any case, the regions where 
the driving stress exceeds the yield stress are expected to 
slide, although the sliding velocity is controlled by mem- 
brane stress connections to nonsliding parts. Because the 
difference between driving and yield stress is large within 
the downstream portions of the weak strips, the sliding will 
initially be very rapid. 

5. Results 

The results in this paper can be reproduced by running a 
script included in the PISM source code distribution. The 
combined time of all of our results is roughly 4600 processor- 
hours on a mix of 2.3 GHz quad core Xeon and 2.6 GHz dual 
core AMD Opteron processors. Runs used 128 processors si- 
multaneously for the 5 km grid results. 
Parameter dependence 

A basic way to illustrate and compare flow is to show the 
surface velocity at a fixed time as in Figure 10. In all cases 
the portions of the weak till strips corresponding to signif- 
icant basal melt produced significant sliding (e.g. > 100 
m/a) at some point in the 5000 year runs, but this snapshot 
at the final time shows some streams are "off". Note that 
in all runs the peak surface speed outside of the weak strips, 
namely the surface speed from the SIA, is approximately 60 
m/a. 

The time evolution of the modeled flow must be illus- 
trated in a different way. A highly- averaged measure is 
the volume of the whole ice sheet. As shown in Figure 11 
the volume decreases in every experiment. The rate of vol- 
ume loss is most rapid in the first 500 model years because 
the high driving stresses in the initial state produce fast 
sliding into the ablation area. This fastest initial flow is 
self-limiting because the geometry changes to produce lower 
driving stresses, and because cold ice is advected down into 
the stream regions, with the ultimate effect that the mod- 
eled basal melt rate decreases or even becomes negative (re- 
freeze), so that the modeled till yield stress eventually in- 
creases at some future time in each stream. 

Experiment PI was continued to a 100 ka run at 10 km 
resolution. The volume time series shown in Figure 12. The 
"leveling out" suggested in the 5 ka result is realized and an 
approximate steady state is achieved. There is sliding flow 
at speeds typical of real ice streams at every time during this 
longer run, as shown in Figures 13 and 14. Those Figures 
show the spatial average of the vertically-integrated horizon- 
tal velocity over the positive thickness, downstream portion 
of each weak strip. There is rapid flow at speeds in ex- 
cess of 10"^ m/a in the first few model years, but these high 
initial speeds are caused by the high driving stresses cre- 
ated by the non-sliding run which created the initial state. 
The modeled flow appropriately re-stabilizes to reasonable 
sliding speeds. The 100 km wide stream in experiment PI 
produces the most volatile flow. The explanation for this is 
that the center of the 100 km strip periodically cools and 



refreezes, causing shutdown of the central part of the sliding 
flow. Oscillations appear in the area-averaged flow speed, 
with period about 900 a (inset in Figure 14). 

The Figures showing flow speeds all include large and 
very brief excursions to much faster flow; these are vertical 
lines sometimes going off scale. These correspond to discrete 
advance and retreat of the margin by a single or a few grid 
spaces. We believe that such discrete margin movement re- 
sults in a brief reduction of backpressure, followed by fast 
flow and advance of the margin to again "press" against 
the few essentially non-sliding ice-filled grid spaces at the 
margin and in the ablation area. 

We now show a sequence of Figures (15,16,17,18,19,20) 
from the highest resolution (5 km) runs. 

Figure 15 shows that the two narrower width streams in 
experiment PI reach a level of relatively steady flow at av- 
erage speeds of roughly 100 m/a within the first 1 ka. By 
contrast, the two wider streams enter cycles with average 
speeds varying from under 10 to over 100 m/a. 

Recall that experiment P2 examines how the flow de- 
pends on the grid orientation of weak strips, relative to the 
cartesian grid directions. As shown in Figure 16, although 
the strength and width of the weak strips are the same, 
the resulting speeds are somewhat different. There is rough 
agreement for the first 500 model years. The range of speed 
variation (20 to 200 m/a) roughly agrees for the whole 5000 
years. All three streams enter similar oscillating states, but 
the cycles are out of phase and differ in period (1300 to 
2100 years). Note that because all three streams are mod- 
eled within the same ice sheet there is some "competition" 
for ice near the divide (dome). Nonetheless be believe grid 
orientation has some effect on flow speed, a result which 
needs to be clarified in future work. 

Experiment P3 evaluated the effect of bed slope on the 
modeled flow. We see in Figure 17 that larger slopes gener- 
ate faster flow speeds and greater variability. The period of 
the cycle is only weakly dependent on bed slope, however. 

Experiment P4 evaluated the effect of till friction angle 
in the downstream region of the till. Figure 18 shows that 
smaller till friction angle yields greater volatility in speed, 
with higher highs and lower lows. The variability of the 
modeled stream flow certainly relates to the availability of 
basal water, which modulates till yield stress (subsection 
2.7). Figure 19 shows the basal water thickness W at the 
end of the 5 ka run for P4. It should be compared to the 
P4 surface speed shown in Figure 10 which shows only one 
active ice stream. The only active stream is the one which is 
underlain by the maximal amount of basal water (W = 2 m), 
which is therefore the stream with the lowest till yield stress. 
Though the yield stress is lowest in the active stream, the 
till friction angle is slightly large in the active stream (3°) 
than in the stream with the lowest friction angle (2°). 

Like Figure 10, Figure 20 also shows surface speed for 
experiment P4, but very near the beginning of the run (at 
5 a). We see rapid sliding as the sheet reduces its driving 
stress to match the (recently and suddenly) reduced basal 
resistance which occurs at the start of all our experiements. 
In fact we see that the fastest flow at this time is in excess 
of 3000 m/a, and the weakest till corresponds directly to the 
lowest till friction angle because basal water thickness has 
not yet lost its initial angular symmetry (e.g. as in Figure 
5). 

The surface elevation evolves though the mass continuity 
equation in all runs. Because the sliding is rapid, surface 
evolution occurs rapidly and is mostly stabilized within the 
first 1 ka. Figure 21 illustrates the fact that the high driving 
stresses of the initial state have evolved, by the end of the 
100 ka run, to lower driving stresses. There is roughly con- 
stant surface slope along the sliding portion of the stream; 
compare Figure 2 in [Joughin et al, 2001] which shows the 
surface profile of along approximate fiowlines in the NE 
Greenland ice stream. 
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Effects of grid refinement 

Figure 22 shows that experiment PI behaves stably un- 
der changes in grid spacing when measured by volume. Over 
the full 5000 year evolution there is no obvious convergence 
toward a fine grid result, but the evidence for convergence is 
clear when we consider only the first 1000 years of the run, 
as in Figure 23. In particular, the 7.5 and 5 km grid results 
remain very close for the first 500 years of the run. Grid 
refinement results for experiments P2, P3, and P4 are very 
similar and are not shown. 

The vertical grid was described in the previous section 
and was used in all runs except for one version (10 km hori- 
zontal grid) wherein the number of grid points in the vertical 
was increased by a factor of two. This is labeled "fine vert" 
in the last two Figures. This change in the vertical grid 
makes little difference in these Figures and to essentially all 
results in this paper. 

Grid refinement produces good behavior in terms of the 
spatial distribution of surface speed in modeled ice streams. 
Figure 25 shows the final state of the 70 km wide stream in 
experiment PI, at four horizontal grid spacings. 
Parallelization 

In a parallel application like PISM, ideally a doubling of 
the number of processors would halve the time to complete 
a run. This is never quite achieved in practice for any par- 
allel application because of the time for interprocess com- 
munication and other "overhead". In our case the parallel 
solution of linear systems is the most difficult part to speed 
up through parallelization, but this is handled in a sophisti- 
cated manner by the PETSc library [Balay and eight others, 
2006]. As shown in Figure 24, a perfect rate of speedup is 
not achieved but parallelization is notably effective anyway. 
The change to "asymptotically almost perfect" appearance 
of the speedup in the range from 16 to 128 processors is 
likely an effect of memory parallelism, wherein the part of 
the simulation handled by a processor becomes smaller as 
the number of processors increases, so that it fits in fast 
cache memory. 

6. Discussion 

The experiments above are for a simplified ice sheet with 
a steady and angularly-symmetric climate (surface mass bal- 
ance and temperature). On the other hand, new coupling 
mechanisms interact in a highly dynamic manner, as appro- 
priate to real ice sheet dynamics, so results are not easily 
summarized. Nonetheless we can list some common features 
of the above results: 

(i) All experiments start with a brief period of very fast 
sliding. In the initial state the driving stress is held fully by 
resistance at the base, but once sliding is allowed the driv- 
ing stress is far out of balance. As shown in Figure 20, for 
example, the surface velocity exceeds 3000 m/a at 5 years 
into experiment P4, but it later stabilizes to speeds typical 
of large ice streams with low bed slope (e.g. Figure 10). 

(ii) After an initial period of rapid ice volume loss, caused 
by fast sliding into the ablation area, the ice sheet volume 
stabilizes. Fitting an exponential decay curve to the data 
shown in Figure 12 gives 22 ka as the exponential time for 
volume decay. The estimated steady volume is 2.106 x 10^ 
km*^ in experiment PI, very slightly smaller than the mean 
steady volume of 2.128 x 10^ km^ in EISMINT II experiment 
A [Payne et al, 2000]. 

(Hi) Surface velocity within ice streams sometimes enters 
a limit cycle. Wider ice streams are more likely to continue 
cycling because the center of the stream is the location at 
which advected cold ice can cause basal refeeze (followed by 
an increase in till yield strength and sliding shut down). In 
a wider stream there is more space for advection to domi- 
nate the energy balance, away from the intense dissipation 
heating at stream margins. 



(iv) Experiments without a bedrock thermal layer (not 
shown) suggest that having a bedrock layer of modest thick- 
ness (e.g. > 100 m) is important. Without it the basal lu- 
brication dynamics are more erratic. The need for a modest 
layer identified here, to mollify short timescale variations 
in basal lubrication, is distinct from the known importance 
of a relatively deep bedrock thermal layer in order to cor- 
rectly assimilate long time scale climate variations into the 
ice sheet temperature field. 

{v) Schoof [2006a] observes that for the plastic till ver- 
sion of the SSA, the stability of time-dependent geometry 
evolution is unknown as a mathematical matter. In fact 
there is no proof that cliffs will not appear on the ice sur- 
face at sliding/nonsliding transition locations in the Schoof 
[2006a] model, if it were to be used by itself in a time- 
dependent manner. Our model adds horizontal plane shear- 
ing by the SI A, which causes the equation of mass continuity 
to be partly diffusive, which smooths the surface. It there- 
fore exhibits stable geometry evolution. 

The fundamental explanation for the variability of the 
stream flow in our model is exactly the mechanism described 
in [Raymond, 2000], in which strain dissipation and sliding 
friction heating compete with a combination of cold ice ad- 
vection and basal drainage, in providing lubrication for slid- 
ing flow. In our case, basal drainage is modeled only through 
limiting the effective basal water thickness to 2 m. 

We note three caveats about our model, beyond the shal- 
lowness assumptions addressed in the introduction: 

(a) Experiments in this paper do not include floating ice 
shelves, which are subject to zero till yield stress and to 
the floatation criterion. On the other hand, the same PISM 
code for solving the SSA has successfully duplicated the EIS- 
MINT diagnostic comparison of modeled to observed Ross 
ice shelf velocities [Bueler et al, 2008; MacAyeal, 1989]. 

(b) Modeling the flow of ice in the immediate vicinity 
of the grounded margins apparently requires more complete 
stress balances than used here. 

(c) Detailed modeling of ice stream shear margins is likely 
to require more complete stress balances than used here. 
Our model allows transmission of side resistance into the 
interior of the sliding ice stream through membrane stresses 
only. Real shear margins involve bridging stresses, exotic 
concentration of strain dissipation heating and crevasse- 
related cooling [Truffer and Echelmeyer , 2003], and accumu- 
lation of anisotropy and damage in the ice fabric, all missing 
in our model. 

These caveats each suggest the importance of high resolu- 
tion, thermomechanically coupled "higher-order" and (full) 
Stokes stress balance modeling, either at whole ice sheet 
scale or in locally-refined submodels within a shallower 
whole sheet model like the one here. 

7. Conclusion 

The open source model in this paper includes membrane 
stress balance, is thermo-mechanically coupled, includes a 
basal strength parameterization depending on basal melt, 
and is high resolution. We demonstrate runs for a Green- 
land scale ice sheet on 5 km grid for a duration of 5000 model 
years, and for 100 ka on a 10 km grid, so a model with these 
features is now available for most prognostic modeling pur- 
poses. 

In fact, our model is of intermediate computational ex- 
pense between the easier thermomechanically-coupled SIA 
and more expensive "higher-order" models. Its implemen- 
tation is fully parallel. Specifically, we believe that the nu- 
merical and scientific computing choices made here are an 
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effective step toward parallel implementation at high reso- 
lution of the Blatter [1995] model with a treatment of basal 
sliding following [Schoof, submitted]. 

Our model has distinct advantages over most existing 
whole ice sheet models which include sliding. The basal 
sliding velocity field here is continuous, unlike the sliding 
field generated by the SIA with a classical temperature- 
dependent sliding law. 

The underlying physical model is supported by recent ob- 
servations [Joughin et aL, 2004b; Tulaczyk et aL, 2000] and 
theory [Schoof ^ 2006a] for ice streams. The resulting flow 
speeds are realistic for those parts of real ice sheets where 
bed gradients are modest, including the NE Greenland and 
Siple coast ice streams. Our model is appropriate to the 
shallow southern margins (lobes) of the Laurentide ice sheet, 
which experienced significant streaming. 

We expect that this model can assimilate the informa- 
tion in basal shear stress maps derived by inverse mod- 
elling, though the precise method for deriving till friction 
angle from inversion-derived basal shear stress and veloc- 
ity remains to be addressed. The heuristic concept that 
till is weaker at locations of lower bed elevation, because 
of a marine history for that bed [Huybrechts and de Wolde, 
1999], can be directly implemented, however. This yields 
a promising model of West Antarctic ice streams within a 
whole Antarctic ice sheet model (in preparation). Short 
time-scale, local changes basal lubrication of the type as- 
sociated to moulin drainage of surface melt [Zwally et a/., 
2002] can also be simulated directly without inverse meth- 
ods, by forcing changes in modeled pore water pressure. 

Appendix A: Finite differences for SSA 
stress balance 

We use the notation of subsection (2.6), but let x = xi 
and y = X2 and we denote an approximation of the function 
value f{x,y) at a grid point [x^^y-^) by f^'-^ . We use differ- 
ence notation (5+,/'^' = f'^^'-p', 6-.P' = P'-f''^', 
and Ax/*'-^ = j^+i.i _ p ^ q^j^^j corresponding notation 
for y differences. Temporarily denoting the product v H hj 
'W, our approximation of the first of Equations (17) is 



Ax 



Ax 



Ax 



4Ay 
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Ay 



4Ax 
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The second of Equations (17) is approximated by a similar 
scheme. In computing the staggered value oi v H the thick- 
ness H is averaged onto the staggered grid, while staggered 
grid values of the effective viscosity T> come from Equation 
(14) using the same centered finite differences as in Equation 
(Al). Figure 26 shows which of the quantities v i7, vi^V2 are 
evaluated at which grid points in our scheme. 

To actually solve finite difference Equation (Al) requires 
an "outer" nonlinear iteration and an "inner" linear iter- 
ation. The outer iteration produces a sequence of sparse 
linear problems in which the values of i>H are (temporarily) 
known. The inner iteration solves the linear problem for vi 
and V2. 

The outer iteration continues until successive values of 
VH are within a specified tolerance. The criterion used for 
the results here is that the norm of the difference of suc- 
cessive values of VH is less than 10""^ of the norm of 



I'M itself. We do not require convergence of the velocities 
themselves in the outer iteration. 

The inner linear problems can each be thought of as 
"Av = b", with the approximated driving stress pgHVh 
as b. The membrane stresses and the basal stress Th all ap- 
pear on the left side as parts of the matrix A. We factor 
Th into a coefficient, which depends in a strongly nonlinear 
way on the velocity, times the velocity. 

If the grid has Mi x M2 horizontal points then there are 
2M1M2 scalar unknowns. Using Mi = M2 = 301 as in the 5 
km grid results shown, each linear system has 1.8 x 10^ rows. 
Five grid values of I'l and eight grid values of V2 are involved 
in Equation (Al), so there are 13 nonzero coefficients per 
row. These sparse linear problems are solved using a par- 
allel iterative solver within PETSc [Balay and eight others, 
2006]. For the results here we used the PETSc default, GM- 
RES(30) with block Jacobi and ILU preconditioning [Saad, 
2003]. This solver /preconditioner pair is just one choice of 
many "Krylov" iterative methods within PETSc, and other 
combinations work. We use the default Krylov relative tol- 
erance of 10~^. 



Appendix B: Traditional SIA sliding laws 
and their difficulties 

Consider a slab on a slope [Paterson, 1994] of thickness 
Ho and slope a, with horizontal coordinate x, positive- 
upward vertical coordinate z, and horizontal velocity u. 

Sliding laws of the type used in EISMINT II experiment 
H [Payne et al, 2000] and ISMIP-HEINO [Greve et al, 2006] 
give the sliding velocity as a function of the driving stress 
when the basal ice reaches the pressure melting tempera- 
ture Tq . Suppose that the basal temperature T increases 
along the fiow, reaching pressure melting temperature Tq at 
some location xo, as shown in Figure 27. Suppose T < Tq 
for X < xq while T = Tq for x > xq. Concretely, suppose 
Ub = when T < Tp and Ub = CpgHQa when T = Tq , for a 
positive constant C. The horizontal velocity at any location 
(x, z) in the ice is 



u{x, z) = Ub{x) + 



/■ 

Jb{x) 



du 



The SIA computes du/dz from the ice temperature, surface 
slope, and depth below the surface. It follows that there is a 
jump in the basal velocity which implies a jump in horizon- 
tal velocity at every vertical level z between the base and 
the surface of the ice: 

u{xq , z) — u(xq ,z) = CpgHQa. 

On the other hand, the ice is incompressible. Therefore 
the vertical velocity is formally infinite at all points in this 
ice column. Indeed, the vertical velocity is 



w(x, z) ■ 



f 



du 
dx 



dC, 



at any place where the horizontal velocity is differentiable 
in X. In the ice column ai x = xq we necessarily interpret 
''du/dx" as infinity. (We have made harmless small-slope 
approximations, including assuming it; = at the base of 
the ice.) 

It remains to explain why this problem has not already 
stopped this ice sheet modeling practice. Consider a numer- 
ical scheme in which xq is between a pair of gridpoints x- 
and x+ which are separated by Ax. Suppose Hq — 2000 
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m, a = 0.001, p = 910 kgm"^ g 9.81 ms"^ and 
C = 10"^mPa"^a"^ [c.f. Payne et al, 2000]. If 

du , . uix-L, z) — u(x-, z) 
ox Ax 

then Table 4 gives estimates of the surface value of verti- 
cal velocity, starting with common values of Ax used in ice 
sheet modeling, and then down to 1 km. We see that for 
the rough grids used in EISMINT and ISMIP-HEINO, for 
example, this "infinity" is well-hidden unless the modeler 
carefully examines locations of slightly anomolous vertical 
velocities. The problem starts appearing at the level of grid 
refinement which is necessary to resolve the stresses in indi- 
vidual ice streams, however. 

Because the vertical velocity is a part of the energy con- 
servation Equation (7), however, anomolous vertical veloc- 
ity is merely the start of unreasonable, time-dependent, 
thermomechanically-coupled behavior. We hypothesize that 
the strange spokes in EISMINT II experiment H [Payne 
et al.^ 2000] are a consequence of the mechanism described 
in this appendix. By contrast, the spokes for the nonslid- 
ing experiments in EISMINT II came from a different ther- 
momechanical fluid instability mechanism, as discussed in 
[Bueler et aL, 2007] and references. 

We also hypothesize that the "Heinrich events" demon- 
strated in ISMIP-HEINO [Greve et al, 2006] are conse- 
quences of thermomechanical coupling to this anomolous 
mechanism here. This comment does not discount the 
possibility that ice sheet sliding is connected to Heinrich 
events. Rather, we believe that continuum models for Hein- 
rich events should not use SIA sliding laws of the type de- 
scribed in this Appendix. Indeed, membrane stresses mod- 
ulate sliding so the need to be included in any model for 
time-dependent ice stream behavior. 
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Table 1. Notation, units, and values for constants. Vectors are in bold. 

Symbol Description SI Units Value 

seconds per year 31556926 

V, V- gradient and divergence in horizontal coordinates 

A{T*) temperature-dependent rate factor in flow law; Eqn (5) Pa~^s~^ 

^(T*) ice hardness; 5(r*) = A(r*)-^/^ Pa s^/^ 

b bed elevation m 

P dependence of melting point on depth Km~^ 8.66 x 10~^ 

Ci specific heat capacity for ice J (kg K)~^ 2009 

Dij strain rate tensor s""*^ 

Dl depth of lithosphere thermal model m 515 

DsiA diffusivity associated to SIA; Eqn (24) m^s~^ 

S regularizing velocity for plastic till; Eqn (27) ms~^ 0.01 ma""*^ 

e regularizing velocity for effective viscosity; Eqn (26) ms~^ 1.0 ma~^ 

/(|v|) weighting function for combining velocities; Eqns (21), (22) 

till friction angle; Eqn (19) 

Go geothermal flux Wm"^ .042 

g acceleration of gravity ms~^ 9.81 

h ice surface elevation m 

H ice thickness m 

ki thermal conductivity of ice W(K m)~^ 2.10 

Kr thermal diffusivity of lithosphere m^s~^ 9.09x10"^ 

L latent heat of fusion for ice Jkg~^ 3.35 x 10^ 

L diffusion distance for melt water thickness; Eqn (11) m 20 km 

Ljy regularizing length for effective viscosity; Eqn (26) m 10^ 

M ice-equivalent surface mass balance (M > is accum.) ms~^ 

n Glen exponent 3 

u vertically- averaged (effective) viscosity in SSA; Eqn (14) Pa s 

Q horizontal ice flux m^s~^ 

p pressure (in ice) Pa 

Pw pore water pressure in till; Eqn (20) Pa 

R gas constant J(molK)~^ 8.314 

Pi density of ice kgm~^ 910 

S ice-equivalent basal mass balance (5* > is melting) ms~^ 

SIA "shallow ice approximation", esp. Eqn (13) 

SSA "shallow shelf approximation", esp. Eqn (17) 

T absolute ice (or bedrock) temperature K 

T* pressure corrected ice temperature; = To — ^{h — z) K 

To melting temperature for ice K 273.15 

Ts surface temperature K 

t time s 

t diffusion time for melt water thickness; Eqn (11) s 1000 a 

Th^i basal shear stress applied to ice; Eqn (18) Pa 

Tc till yield stress Pa 

Tij deviatoric stress tensor Pa 

u = {ui^U2^uz) velocity computed from non-sliding SIA; Eqn (13) ms~^ 

U = (t/i, U2, Us) ice velocity; note Eqn (21) ms~^ 

Jj vertically-averaged horizontal velocity ms~^ 

JJb basal value of the ice velocity ms~^ 

V = {vi^V2,V3) velocity computed from plastic till SSA; Eqn (17) ms~^ 

W effective thickness of stored basal water; Eqn (11) m 

Wo limit on basal drainage parameter m 2.0 

xi,X2,X3 cartesian coordinates m 

z alternate notation for X3; positive upward m 

Table 2. Constants for the exact SSA solution which is equation 4.3 in [Schoof, 2006a], with values B,f from the same source. 

Symbol Meaning Value 

B ice hardness 3.7 x lO^Pas^/^ 

/ scale for till yield stress 17.854 kPa 

ho constant ice thickness 2000 m 

L half-width of weak till 40 km 

m power 10 

tan6> bed slope 0.001 
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Table 3. Parameter studies: Each experiment runs for 5000 
models years on grids of 15, 10, 7.5, and 5 km spacing, starting 
from the same initial state. 

Name Parameter explored Comments 

PI width of weak till strip 10 km grid run extended to 100k years 

P2 orientation of weak till strip only three weak till strips 

P3 bed slope 

P4 strength of downstream till 



Table 4. Numerical surface value of the vertical velocity from a traditional SIA sliding law. 
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Figure 1. A hierarchy of ice dynamics models (ellipses) 
with sliding parameterizations which have been applieds 
to the shallower models (boxes). Solid arrows show rig- 
orous small-parameter shallowness arguments, while Uh 
denote basal ice velocity and Th basal shear stress. 
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Figure 3. Solid curve: Ice stream velocity v{y)m the 
down slope direction in the exact SSA solution. Dot- 
ted: Till yield stress Tc{y) grows sharply at \y\ =40 km, 
and so significant sliding occurs only within the inter- 
val —40 km < y < 40 km. The region of sliding is not 
predetermined. 
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Figure 4. Numerical errors in along-flow velocity for 
exact SSA solution shown in Figure 3. Convergence at 
rate 0(Ay^'^°) is found by fitting errors for grids with 
160 <Ay< 5000 m. 
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Figure 5. Ice thickness (left) and thickness (right) of 
basal water layer for the initial state. Effective basal 
water thickness W ranges from zero (blue) to 2 m (red); 
the central blue has frozen base and zero basal water. 
The right figure, and all other color map-plane Figures 
in this section, show a 1500 x 1500 km square region. 




Figure 6. Maps of till friction angle (j){x,y) for exper- 
iments PI (a), P2 (b), and P4 (c). Red is 15°. Blue in 
weak strips in a,b is 5°, also in upstream parts of weak 
strips in c. See text for downstream values in experi- 
ment P4 and for description of experiment P3. 



20 




30 (km) 

--50 (km) 
70 (km) 



■■■ 100 (km) 

%0 -60 -40 -20 20 40 60 80 
distance from center of strip (km) 

Figure 7. Till friction angle (j) across the four weak till strips in experiment PI, with given strip widths. 
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Figure 8. Trough bed topography in experiment P3. 
Red is flat plateau, and troughs descend m (east), 500 
m (north), 1000 m (west), and 2000 m (south) below the 
plateau. 



Figure 9. Difference At = pgH\S/h\ — Tc between basal 
driving stress and till yield stress at start of experiment 
PI. Tc computed from Equation (19). Darkest red re- 
gions have At > 5 x 10^ Pa, so they slide strongly at the 
beginning of the run. 
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Figure 10. Modeled horizontal surface ice speed (m/a) 
in logarithmic color scale, at 5 km resolution and at the 
end of 5 ka runs: experiments PI (upper left), P2 (upper 
right), P3 (lower left), P4 (lower right). 
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Figure 11. Evolution of model ice sheet volume for the four experiments at 5 km grid resolution. 
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Figure 12. Evolution of model ice sheet volume for experiment PI at 10 km grid resolution over a 100 ka span. 
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Figure 13. Average downstream speed (magnitude of 
vertically-averaged horizontal velocity) for three of the 
four ice streams in experiment PI, over a 100 ka span. 
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Figure 14. Average downstream speed for the 100 km 
wide ice stream in experiment PI, over a 100 ka span. 
Inset shows final 5 ka of run. 
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Figure 15. Average downstream speed for the ice 
streams in experiment PI, at 5 km resolution for a 5 
ka run. Weak till strip width is given in the legend. 
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Figure 16. Average downstream speed for the three 
ice streams in experiment P2, at 5 km resolution. Strip 
angle, relative to the closest grid (cardinal) direction, is 
given in the legend. 
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Figure 17. Average downstream speed for the ice 
streams in experiment P3, at 5 km resolution. Trough 
slope is given in the legend as a percentage. 



2^' 




t (a) 



Figure 18. Average downstream speed for the ice 
streams in experiment P4, at 5 km resolution. Down- 
stream till friction angle is given in the legend. 




Figure 19. Modeled basal water effective thickness (m) 
for experiment P4, at 5 km resolution and at the end 
of the 5 ka run. The effective thickness is limited to a 
maximum of 2 m in the model. 
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Figure 20. Modeled horizontal surface ice speed (m/a) 
from experiment P4, at 5 km resolution but at only 5 
model years from the start of the run. 
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Figure 21. Surface elevation, from the center of the ice 
sheet along the flow line corresponding to the center of 
the 70 km wide weak till strip in experiment PI, at 10 
km resolution and at the end of the 100 ka run, compared 
to the initial state. 
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Figure 22. Effect of grid refinement on the volume time series for the 5 ka run of experiment PI. 
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Figure 23. Effect of grid refinement on the volume time series in first 1 ka (experiment PI). 
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Figure 25. Detail of modeled horizontal surface ice 
speed (m/a) for experiment PI, in the same color scale 
used for Figure 10, at the end of 5 ka runs. Horizontal 
resolutions (top to bottom): 15, 10, 7.5, and 5 km. 
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Figure 26. A stencil for finite difference scheme (Al), 
for approximating the first of Equations (17). Triangles 
show staggered grid points where i>H is evaluated. Cir- 
cles and squares show where vi and V2 are approximated, 
respectively. 
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Figure 27. If the sliding law "turns on" at a point where 
the pressure-melting temperature T — Tq then there is 
a discontinuity in the sliding velocity everywhere in the 
ice column above that point. The grey band suggests 
horizontal velocity it as a function of elevation. 



